#include "OclLapack.h"
#include <cmath>
#include <algorithm>
#include <vector>
#include <memory.h>
#include "../BindKernelArgs.h"
#include "../ClMem.h"
#include <limits>

using namespace std;

#define COUNT(x) (sizeof(x) / sizeof(x[0]))

extern "C" {
	int lsame_(char *, char *);
	int xerbla_(char *, int *);
}

extern int ilaenv_(int *, const char *, const char *, int *, int *, int *, int *);
extern double slamch_(char *);

enum Range {
	RANGE_NONE = 0,
	RANGE_ALL,
	RANGE_VALUE,
	RANGE_INDEX,
};

enum Order {
	ORDER_NONE = 0,
	ORDER_ENTIRE_MATRIX,
	ORDER_BY_BLOCK,
};

struct Interval {
	float low;
	float up;
};

Interval getGershgorinInterval(int n, float *d, float *e) {
	Interval g;
	g.low = g.up = d[0];
		
	float prev = 0;
	for (int j = 0; j < n - 1; ++j) {
		float absE = fabs(e[j]);
		g.low = min(g.low, d[j] - prev - absE);
		g.up = max(g.up, d[j] + prev + absE);		
		prev = absE;
	}

	g.low = min(g.low, d[n - 1] - prev);
	g.up = max(g.up, d[n - 1] + prev);	
	return g;
}

template <typename T>
inline T sqr(T x) {
	return x * x;
}

// Compute N(c), the number of eigenvalues less than c
int getNumEigensLess(float c, int n, float *d, float *e, float pivmin) {
	
	float diff = d[0] - c;
	diff = select(diff, -pivmin, fabs(diff) < pivmin);
	int count = select(0, 1, diff <= 0);

	for (int j = 1; j < n; ++j) {
		diff = d[j] - c - sqr(e[j - 1]) / diff;
		diff  = select(diff, -pivmin, fabs(diff) < pivmin);
		count += select(0, 1, diff <= 0);
	}
	return count;
}

int isNarrow(Interval interval, float abstol, float reltol) {
	float width = fabs(interval.up - interval.low);
	float norm = max(fabs(interval.up), fabs(interval.low));
	float minWidth = max(abstol, reltol * norm);
	return (width < minWidth);
}

ClEvent OclLapack::recalculateEigenIntervals(size_t intervalCount, int n,
		cl_mem d, cl_mem e,		
		cl_mem inIntervals,	cl_mem outIntervals,
		cl_mem inNumLess, cl_mem outNumLess,		
		float abstol, float reltol, float pivmin) 
{
	BindKernelArgs(recalculateEigenInterval) << n <<
			d << e <<
			inIntervals << outIntervals <<
			inNumLess << outNumLess <<
			abstol << reltol << pivmin;
	
	cl_event event;
	CL_CHECK(clEnqueueNDRangeKernel(queue,
			recalculateEigenInterval, 1, NULL,
			&intervalCount, NULL,
			0, NULL, &event));
	return ClEvent(event);
}

bool isAllNarrow(int count, Interval *intervals, float abstol, float reltol) {
	for (int j = 0; j < count; j++) {
		if (!isNarrow(intervals[j], abstol, reltol)) return false;		
	}
	return true;
}

ClEvent OclLapack::calculateNumEigensInIntervals(size_t numEigens,
	cl_mem numEigensInInterval, cl_mem intervals,
	int n, cl_mem d, cl_mem e,
	float pivmin) 
{
	BindKernelArgs(calculateNumEigensInInterval)
			<< numEigensInInterval << intervals
			<< n << d << e
			<< pivmin;
	
	cl_event event;
	CL_CHECK(clEnqueueNDRangeKernel(queue,
			calculateNumEigensInInterval, 1, NULL,
			&numEigens, NULL,
			0, NULL, &event));
	return ClEvent(event);
}

ClEvent OclLapack::recalculateIntervals(size_t numEigens,
		int n, cl_mem d, cl_mem e,
		cl_mem inIntervals, cl_mem outIntervals,
		cl_mem numEigensInIntervals,
		float abstol, float reltol, float pivmin) 
{
	BindKernelArgs(recalculateInterval)
			<< n << d << e 
			<< inIntervals << outIntervals
			<< numEigensInIntervals
			<<	abstol << reltol << pivmin;
	
	cl_event event;
	CL_CHECK(clEnqueueNDRangeKernel(queue,
			recalculateInterval, 1, NULL,
			&numEigens, NULL,
			0, NULL, &event));
	return ClEvent(event);
}

int OclLapack::convergeIntervals(int n, int maxIterations, 
		float *d, float *e,
		void* _intervals,
		int *numLess,
		float abstol, float reltol, float pivmin) 
{	
	Measure measure("convergeIntervals");
	
	Interval *intervals = (Interval*)_intervals;
	
	abstol = max(abstol, pivmin);
	numLess[0] = getNumEigensLess(intervals[0].low, n, d, e, pivmin);
	numLess[1] = getNumEigensLess(intervals[0].up, n, d, e, pivmin);
	int numEigens = numLess[1] - numLess[0];	
	
	ClMem<float> bufferD(queue, n, d);
	ClMem<float> bufferE(queue, n - 1, e);
	ClMem<Interval> currentIntervals(queue, numEigens);
	ClMem<Interval> nextIntervals(queue, numEigens);
	ClMem<int> currentLess(queue, numEigens + 1);
	ClMem<int> nextLess(queue, numEigens + 1);
		
	WAIT_ALL(currentLess.beginWrite(numLess, 0, 2),
		nextLess.beginWrite(numLess, 0, 2),
		currentIntervals.beginWrite(intervals, 0, 1));			
		
	for (int i = 0; i < maxIterations; ++i) {
						
		// recalculate intervals		
		recalculateEigenIntervals(numEigens, n,
				bufferD, bufferE,				
				currentIntervals, nextIntervals,
				currentLess, nextLess,
				abstol, reltol, pivmin).wait();
					
		swap(currentIntervals, nextIntervals);
		swap(currentLess, nextLess);
		
		// read intervals and num less
		currentIntervals.read(intervals, 0, numEigens);
						
		if (isAllNarrow(numEigens, intervals, abstol, reltol)) break;
	}
	
	currentLess.read(&numLess[1], 1, numEigens);
	return numEigens;
}

void findEigenvalue(int n, float *d, float *e, Interval &interval, 
		int index, int maxIterations,
		float abstol, float reltol, float pivmin) 
{	
	int numEigensLessLow = -1;
	int numEigensLessUp = n + 1;
	
	for (int i = 0; i < maxIterations; ++i) {
		if (isNarrow(interval, max(abstol, pivmin), reltol) ||
				 numEigensLessLow >= numEigensLessUp) break;

		float middle = (interval.low + interval.up) * 0.5f;
		int numEigensLessMiddle = getNumEigensLess(middle, n, d, e, pivmin);

		// IJOB=3: Binary search.  Keep only the interval 
		// containing  w  s.t. N(w) = NVAL 
		if (numEigensLessMiddle <= index) {
			interval.low = middle;
			numEigensLessLow = numEigensLessMiddle;
		}
		if (numEigensLessMiddle >= index) {
			interval.up = middle;
			numEigensLessUp = numEigensLessMiddle;
		}
	}
}

int OclLapack::sstebz(char *range, char *order, int *n, float *vl,
		float *vu, int *il, int *iu, float *abstol, float *d, float *e,
		int *m, int *nsplit, float *w, int *iblock, int *
		isplit, float *work, int *iwork, int *info) 
{
	//   -- LAPACK routine (version 3.2) -- 
	//      Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 
	//      November 2006 
	//      8-18-00:  Increase FUDGE factor for T3E (eca) 

	//      .. Scalar Arguments .. 
	//      .. 
	//      .. Array Arguments .. 
	//      .. 

	//   Purpose 
	//   ======= 

	//   SSTEBZ computes the eigenvalues of a symmetric tridiagonal 
	//   matrix T.  The user may ask for all eigenvalues, all eigenvalues 
	//   in the half-open interval (VL, VU], or the IL-th through IU-th 
	//   eigenvalues. 

	//   To avoid overflow, the matrix must be scaled so that its 
	//   largest element is no greater than overflow**(1/2) * 
	//   underflow**(1/4) in absolute value, and for greatest 
	//   accuracy, it should not be much smaller than that. 

	//   See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 
	//   Matrix", Report CS41, Computer Science Dept., Stanford 
	//   University, July 21, 1966. 

	//   Arguments 
	//   ========= 

	//   RANGE   (input) CHARACTER*1 
	//           = 'A': ("All")   all eigenvalues will be found. 
	//           = 'V': ("Value") all eigenvalues in the half-open interval 
	//                            (VL, VU] will be found. 
	//           = 'I': ("Index") the IL-th through IU-th eigenvalues (of the 
	//                            entire matrix) will be found. 

	//   ORDER   (input) CHARACTER*1 
	//           = 'B': ("By Block") the eigenvalues will be grouped by 
	//                               split-off block (see IBLOCK, ISPLIT) and 
	//                               ordered from smallest to largest within 
	//                               the block. 
	//           = 'E': ("Entire matrix") 
	//                               the eigenvalues for the entire matrix 
	//                               will be ordered from smallest to 
	//                               largest. 

	//   N       (input) INTEGER 
	//           The order of the tridiagonal matrix T.  N >= 0. 

	//   VL      (input) REAL 
	//   VU      (input) REAL 
	//           If RANGE='V', the low and up bounds of the interval to 
	//           be searched for eigenvalues.  Eigenvalues less than or equal 
	//           to VL, or greater than VU, will not be returned.  VL < VU. 
	//           Not referenced if RANGE = 'A' or 'I'. 

	//   IL      (input) INTEGER 
	//   IU      (input) INTEGER 
	//           If RANGE='I', the indices (in ascending order) of the 
	//           smallest and largest eigenvalues to be returned. 
	//           1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 
	//           Not referenced if RANGE = 'A' or 'V'. 

	//   ABSTOL  (input) REAL 
	//           The absolute tolerance for the eigenvalues.  An eigenvalue 
	//           (or cluster) is considered to be located if it has been 
	//           determined to lie in an interval whose width is ABSTOL or 
	//           less.  If ABSTOL is less than or equal to zero, then ULP*|T| 
	//           will be used, where |T| means the 1-norm of T. 

	//           Eigenvalues will be computed most accurately when ABSTOL is 
	//           set to twice the underflow threshold 2*SLAMCH('S'), not zero. 

	//   D       (input) REAL array, dimension (N) 
	//           The n diagonal elements of the tridiagonal matrix T. 

	//   E       (input) REAL array, dimension (N-1) 
	//           The (n-1) off-diagonal elements of the tridiagonal matrix T. 

	//   M       (output) INTEGER 
	//           The actual number of eigenvalues found. 0 <= M <= N. 
	//           (See also the description of INFO=2,3.) 

	//   NSPLIT  (output) INTEGER 
	//           The number of diagonal blocks in the matrix T. 
	//           1 <= NSPLIT <= N. 

	//   W       (output) REAL array, dimension (N) 
	//           On exit, the first M elements of W will contain the 
	//           eigenvalues.  (SSTEBZ may use the remaining N-M elements as 
	//           workspace.) 

	//   IBLOCK  (output) INTEGER array, dimension (N) 
	//           At each row/column j where E(j) is zero or small, the 
	//           matrix T is considered to split into a block diagonal 
	//           matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which 
	//           block (from 1 to the number of blocks) the eigenvalue W(i) 
	//           belongs.  (SSTEBZ may use the remaining N-M elements as 
	//           workspace.) 

	//   ISPLIT  (output) INTEGER array, dimension (N) 
	//           The splitting points, at which T breaks up into submatrices. 
	//           The first submatrix consists of rows/columns 1 to ISPLIT(1), 
	//           the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), 
	//           etc., and the NSPLIT-th consists of rows/columns 
	//           ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. 
	//           (Only the first NSPLIT elements will actually be used, but 
	//           since the user cannot know a priori what value NSPLIT will 
	//           have, N words must be reserved for ISPLIT.) 

	//   WORK    (workspace) REAL array, dimension (4*N) 

	//   IWORK   (workspace) INTEGER array, dimension (3*N) 

	//   INFO    (output) INTEGER 
	//           = 0:  successful exit 
	//           < 0:  if INFO = -i, the i-th argument had an illegal value 
	//           > 0:  some or all of the eigenvalues failed to converge or 
	//                 were not computed: 
	//                 =1 or 3: Bisection failed to converge for some 
	//                         eigenvalues; these eigenvalues are flagged by a 
	//                         negative block number.  The effect is that the 
	//                         eigenvalues may not be as accurate as the 
	//                         absolute and relative tolerances.  This is 
	//                         generally caused by unexpectedly inaccurate 
	//                         arithmetic. 
	//                 =2 or 3: RANGE='I' only: Not all of the eigenvalues 
	//                         IL:IU were found. 
	//                         Effect: M < IU+1-IL 
	//                         Cause:  non-monotonic arithmetic, causing the 
	//                                 Sturm sequence to be non-monotonic. 
	//                         Cure:   recalculate, using RANGE='A', and pick 
	//                                 out eigenvalues IL:IU.  In some cases, 
	//                                 increasing the PARAMETER "FUDGE" may 
	//                                 make things work. 
	//                 = 4:    RANGE='I', and the Gershgorin interval 
	//                         initially used was too small.  No eigenvalues 
	//                         were computed. 
	//                         Probable cause: your machine has sloppy 
	//                                         floating-point arithmetic. 
	//                         Cure: Increase the PARAMETER "FUDGE", 
	//                               recompile, and try again. 

	//   Internal Parameters 
	//   =================== 

	//   RELFAC  REAL, default = 2.0e0 
	//           The relative tolerance.  An interval (a,b] lies within 
	//           "relative tolerance" if  b-a < RELFAC*ulp*max(|a|,|b|), 
	//           where "ulp" is the machine precision (distance from 1 to 
	//           the next larger floating point number.) 

	//   FUDGE   REAL, default = 2 
	//           A "fudge factor" to widen the Gershgorin intervals.  Ideally, 
	//           a value of 1 should work, but on machines with sloppy 
	//           arithmetic, this needs to be larger.  The default for 
	//           publicly released versions should be large enough to handle 
	//           the worst machine around.  Note that this has no effect 
	//           on accuracy of the solution. 
	
	Measure measure("sstebz-ocl");
	
	//  Table of constant values 
	static int c__1 = 1;
	static int c_n1 = -1;
	
	//  Local variables 
	int j, ie, je, nb;
	int numEigens, in;
	int iw;
	float wl, wu;
	int nwl;
	float ulp, wlu, wul;
	int nwu;
	int iend, ioff, itmp1, jdisc;
	float atoli;
	int iwoff;	
	int itmax;
	float wkill, rtoli;
	int irange, idiscl;
	float safemn;
	int idiscu;
	int iorder;
	int ncnvrg;
	float pivmin;
	int toofew;

	// Function Body 
	*info = 0;

	// Decode RANGE 
	if (lsame_(range, "A")) {
		irange = RANGE_ALL;
	} else if (lsame_(range, "V")) {
		irange = RANGE_VALUE;
	} else if (lsame_(range, "I")) {
		irange = RANGE_INDEX;
	} else {
		irange = RANGE_NONE;
	}

	// Decode ORDER 
	if (lsame_(order, "B")) {
		iorder = ORDER_BY_BLOCK;
	} else if (lsame_(order, "E")) {
		iorder = ORDER_ENTIRE_MATRIX;
	} else {
		iorder = ORDER_NONE;
	}

	// Check for Errors 

	if (irange <= 0) {
		*info = -1;
	} else if (iorder <= 0) {
		*info = -2;
	} else if (*n < 0) {
		*info = -3;
	} else if (irange == RANGE_VALUE) {
		if (*vl >= *vu) {
			*info = -5;
		}
	} else if (irange == RANGE_INDEX && (*il < 1 || *il > max(1, *n))) {
		*info = -6;
	} else if (irange == RANGE_INDEX && (*iu < min(*n, *il) || *iu > *n)) {
		*info = -7;
	}

	if (*info != 0) {
		int temp = -(*info);
		xerbla_("SSTEBZ", &temp);
		return 0;
	}

	// Initialize error flags 

	*info = 0;
	ncnvrg = false;
	toofew = false;

	// Quick return if possible 

	*m = 0;
	if (*n == 0) return 0;

	// Simplifications: 

	if (irange == RANGE_INDEX && *il == 1 && *iu == *n) {
		irange = RANGE_ALL;
	}

	// Get machine constants 
	// NB is the minimum vector length for vector bisection, or 0 
	// if only scalar is to be done. 

	safemn = 1 / numeric_limits<float>::max();
	ulp = slamch_("P");
	rtoli = ulp * 2.f;
	nb = ilaenv_(&c__1, "SSTEBZ", " ", n, &c_n1, &c_n1, &c_n1);
	if (nb <= 1) {
		nb = 0;
	}

	// Special Case when N=1 

	if (*n == 1) {
		*nsplit = 1;
		isplit[0] = 1;
		if (irange == RANGE_VALUE && (*vl >= d[0] || *vu < d[0])) {
			*m = 0;
		} else {
			w[0] = d[0];
			iblock[0] = 1;
			*m = 1;
		}
		return 0;
	}

	// Compute Splitting Points 

	*nsplit = 1;
	work[-1 + *n] = 0.f;
	pivmin = 1;

	for (int j = 1; j < *n; ++j) {
		float sqrE = sqr(e[j - 1]);
		if (fabs(sqr(d[j] * ulp)) + safemn > sqrE) {
			isplit[*nsplit - 1] = j;
			++(*nsplit);
			work[j - 1] = 0;
		} else {
			work[j - 1] = sqrE;
			pivmin = max(pivmin, sqrE);
		}
	}
	isplit[*nsplit - 1] = *n;
	pivmin *= safemn;

	// Compute Interval and ATOLI	
	if (irange == RANGE_INDEX) {

		// RANGE='I': Compute the interval containing eigenvalues 
		// IL through IU. 

		// Compute Gershgorin interval for entire (split) matrix 
		// and use it as the initial interval
		Interval g = getGershgorinInterval(*n, d, e);
		float norm = max(fabs(g.low), fabs(g.up));
		g.low -= norm * 2.1f * ulp * *n + pivmin * 4.2000000000000002f;
		g.up += norm * 2.1f * ulp * *n + pivmin * 2.1f;

		// Compute Iteration parameters 
		itmax = (int) ((log(norm + pivmin) - log(pivmin)) / log(2.f)) + 2;
		atoli = (*abstol <= 0) ? ulp * norm : *abstol;

		Interval left = g;
		findEigenvalue(*n, d, e, left, *il - 1, itmax, atoli, rtoli, pivmin);
		
		Interval right = g;
		findEigenvalue(*n, d, e, right, *iu, itmax, atoli, rtoli, pivmin);
				
		wl = left.low;
		wlu = left.up;

		wul = right.low;			
		wu = right.up;		
	} else {
		float norm = max(fabs(d[0]) + fabs(e[0]), fabs(d[*n - 1]) + fabs(e[*n - 2]));
		for (int j = 1; j < *n - 1; ++j) {
			norm = max(norm, fabs(d[j]) + fabs(e[j - 1]) + fabs(e[j]));
		}
		
		atoli = max(*abstol, ulp * norm);		

		if (irange == RANGE_VALUE) {
			wl = *vl;
			wu = *vu;
		} else {
			wl = 0;
			wu = 0;
		}
	}

	// Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. 
	// NWL accumulates the number of eigenvalues .le. WL, 
	// NWU accumulates the number of eigenvalues .le. WU 

	*m = 0;
	iend = 0;
	*info = 0;
	nwl = 0;
	nwu = 0;

	for (int block = 1; block <= *nsplit; ++block) {
		ioff = iend;
		int ibegin = ioff;
		iend = isplit[-1 + block];
		in = iend - ioff;

		if (in == 1) {

			// Special Case -- IN=1 

			if (irange == RANGE_ALL || wl >= d[ibegin] - pivmin) {
				++nwl;
			}
			if (irange == RANGE_ALL || wu >= d[ibegin] - pivmin) {
				++nwu;
			}
			if (irange == RANGE_ALL || wl < d[ibegin] - pivmin && wu >= d[ibegin] - pivmin) {
				++(*m);
				w[-1 + *m] = d[ibegin];
				iblock[-1 + *m] = block;
			}
		} else {

			// General Case -- IN > 1 

			// Compute Gershgorin Interval and use it as the initial interval 
			Interval g = getGershgorinInterval(iend - ibegin - 1, &d[ibegin], &e[ibegin]);
			float norm = max(fabs(g.low), fabs(g.up));
			g.low -= norm * 2.1f * ulp * in + pivmin * 2.1f;
			g.up += norm * 2.1f * ulp * in + pivmin * 2.1f;

			// Compute ATOLI for the current submatrix 
			atoli = max(ulp * max(fabs(g.low), fabs(g.up)), *abstol);

			if (irange == RANGE_INDEX || irange == RANGE_VALUE) {
				if (g.up < wl) {
					nwl += in;
					nwu += in;
					continue;
				}
				g.low = max(g.low, wl);
				g.up = min(g.up, wu);
				if (g.low >= g.up) {
					continue;
				}
			}

			Interval *intervals = (Interval*)&work[*n];
			int *numEigensLess = &iwork[0];
			
			int itmax = (int) ((log(g.up - g.low + pivmin) - log(pivmin)) / log(2.f)) + 2;
			intervals[0] = g;
						
			int numEigens = convergeIntervals(in, itmax, 
					&d[ibegin], &e[ibegin],
					&intervals[0], &numEigensLess[0],
					atoli, rtoli, pivmin); 
									
			nwl += numEigensLess[0];
			nwu += numEigensLess[numEigens];
			iwoff = *m - numEigensLess[0];

			// fill eigenvalues and convergence flags
			for (int i = 0; i < numEigens; i++) {
				float eigen = (intervals[i].low + intervals[i].up) * 0.5f;
			
				bool converged = isNarrow(intervals[i], max(atoli, pivmin), rtoli);
				ncnvrg = ncnvrg || (!converged);
				int ib = converged ? block : -block;
				
				int start = numEigensLess[i] + iwoff;
				int count = numEigensLess[i + 1] + iwoff;				
				for (int j = start; j < count; ++j) {
					w[j] = eigen;
					iblock[j] = ib;
				}
			}			

			*m += numEigens;
		}
	}

	// If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU 
	// If NWL+1 < IL or NWU > IU, discard extra eigenvalues. 

	if (irange == RANGE_INDEX) {
		numEigens = 0;
		idiscl = *il - 1 - nwl;
		idiscu = nwu - *iu;

		if (idiscl > 0 || idiscu > 0) {
			for (je = 1; je <= *m; ++je) {
				if (w[-1 + je] <= wlu && idiscl > 0) {
					--idiscl;
				} else if (w[-1 + je] >= wul && idiscu > 0) {
					--idiscu;
				} else {
					++numEigens;
					w[-1 + numEigens] = w[-1 + je];
					iblock[-1 + numEigens] = iblock[-1 + je];
				}
			}
			*m = numEigens;
		}
		if (idiscl > 0 || idiscu > 0) {

			// Code to deal with effects of bad arithmetic: 
			// Some low eigenvalues to be discarded are not in (WL,WLU], 
			// or high eigenvalues to be discarded are not in (WUL,WU] 
			// so just kill off the smallest IDISCL/largest IDISCU 
			// eigenvalues, by simply finding the smallest/largest 
			// eigenvalue(s). 

			// (If N(w) is monotone non-decreasing, this should never 
			// happen.) 

			if (idiscl > 0) {
				wkill = wu;
				for (jdisc = 1; jdisc <= idiscl; ++jdisc) {
					iw = 0;
					for (je = 1; je <= *m; ++je) {
						if (iblock[-1 + je] != 0 && (w[-1 + je] < wkill || iw == 0)) {
							iw = je;
							wkill = w[-1 + je];
						}
					}
					iblock[-1 + iw] = 0;
				}
			}
			if (idiscu > 0) {

				wkill = wl;
				for (jdisc = 1; jdisc <= idiscu; ++jdisc) {
					iw = 0;
					for (je = 1; je <= *m; ++je) {
						if (iblock[-1 + je] != 0 && (w[-1 + je] > wkill || iw == 0)) {
							iw = je;
							wkill = w[-1 + je];
						}
					}
					iblock[-1 + iw] = 0;
				}
			}
			numEigens = 0;
			for (je = 1; je <= *m; ++je) {
				if (iblock[-1 + je] != 0) {
					++numEigens;
					w[-1 + numEigens] = w[-1 + je];
					iblock[-1 + numEigens] = iblock[-1 + je];
				}
			}
			*m = numEigens;
		}
		if (idiscl < 0 || idiscu < 0) {
			toofew = true;
		}
	}

	// If ORDER='B', do nothing -- the eigenvalues are already sorted 
	// by block. 
	// If ORDER='E', sort the eigenvalues from smallest to largest 

	if (iorder == ORDER_ENTIRE_MATRIX && *nsplit > 1) {
		for (je = 1; je < *m; ++je) {
			ie = 0;
			float sqrE = w[-1 + je];
			for (j = je + 1; j <= *m; ++j) {
				if (w[-1 + j] < sqrE) {
					ie = j;
					sqrE = w[-1 + j];
				}
			}

			if (ie != 0) {
				itmp1 = iblock[-1 + ie];
				w[-1 + ie] = w[-1 + je];
				iblock[-1 + ie] = iblock[-1 + je];
				w[-1 + je] = sqrE;
				iblock[-1 + je] = itmp1;
			}
		}
	}

	*info = 0;
	if (ncnvrg) {
		++(*info);
	}
	if (toofew) {
		*info += 2;
	}
	return 0;

	// End of SSTEBZ 

} // sstebz_ 